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ABSTRACT 


A new approach to numerically solving the problem of 
the constricted electric arcjet is presented. An Euler 
Implicit finite difference scheme is used to solve the full 
compressible Navier Stokes equations in two dimensions. 

The boundary and initial conditions represent the constric- 
tor section of the arcjet and hydrogen is used as a propel- 
lant. The arc is modeled as a gaussian distribution across 
the centerline of the constrictor. 

Temperature, pressure and velocity profiles for 
steady state converged solutions show both axial and radial 
changes in distributions resulting from their interaction 
with the arc energy source for specific input conditions. 
The temperature rise is largest at the centerline where 
there is the greatest concentration of arc energy. The 
solution does not converge for all initial inputs and the 
limitations in the range of obtainable solutions are 
discussed. 
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Chapter I 


INTRODUCTION 

1.1 Review of Arcjet Propulsion 

The concept of electric propulsion has been intro- 
duced as an alternative to the chemical liquid and solid 
rocket. The idea was an outgrowth of the increasingly 
popular arc heaters being used in the late fifties and 
early sixties in wind tunnels for reentry simulation. One 
of the most significant advantages to using electric pro- 
pulsion is the tremendous increase in specific impulse 
(thrust per unit mass flow) over conventional rockets. On 
the average, the electric thruster can get ten times more 
specific impulse than the solid rocket engine. The main 
focus of arcjet research is for orbit -maneuvering applica- 
tions. The high specific impulse means less fuel has to be 
boosted into orbit and that high payload capabilities are 
promising. The arcjet can generate low levels of thrust 
for long periods of time, but in general the transfer times 
for these systems tend to be longer than for other propul- 
sive systems.^ 

Various configurations of electric arc heaters 
exist, but the one more frequently used for propulsive 


\ 


1 




2 


purposes is the concept of the constricted arc (see Fig. 

1). The arc itself is an electric discharge capable of 
sustaining large currents between two electrodes, from the 
negative cathode to the positive anode. The arc is con- 
tained in a narrow constricted duct which keeps it from 
dissipating before it reaches the diverging nozzle. The 
purpose of the constrictor is to provide stability and 

restrict the movement of the arc to prevent radial kink- 
2 

ing. Gas enters axially into the constrictor from either 
side of the cathode, absorbing energy from the arc and 
accelerating as it flows through the section. The moving 
gas keeps the arc from attaching to the walls . Negligible 
heat conduction takes place between the edges of the arc 
and the walls, therefore the cool layer of gas lets the 
wall temperatures remain fixed. The gas expands in the 
divergent nozzle at the constrictor exit thereby producing 
thrust. The arc attaches to the anode at the downstream 
end of the constrictor. 

Previous work with propellants reveals that the 
lower the molecular weight, the higher the exhaust velocity 
that can be obtained.'*' A propellant with low molecular 
weight produces higher levels of specific impulse, such as 
those being used in the current analysis. At these levels, 
the maximum arc chamber temperature using hydrogen will be 
lower than for any other known propellant to obtain the 
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same specific impulse. Although still not altogether 
practical because of storage problems, hydrogen is consid- 
ered an efficient and desirable propellant for applications 
of this kind. In addition, accurate data of hydrogen gas 
properties are readily available and are not considered to 
be a significant source of error in numerical calculations. 

1.2 Previous Theoretical Work 

In the early sixties there was a considerable effort 
made to solve, both experimentally and theoretically, the 
problem of the wall-constricted, direct current electric 
arc. Several analytical models were developed before 
numerical solutions were available. Although the linear- 
ization of the equations and the necessary assumptions that 
must be made for the analytical solution limit their 
effectiveness, significant overall trends can be determined 
and are helpful as a first approximation. 

The simplest analytical solution to the constricted 

4 

arc problem is the Stine-Watson model. In 1962, H. Stine 
and V. Watson modeled an arc air-heater, linearizing and 
uncoupling the energy equation so that it could be solved 
using modified Bessel Functions in one dimension. The 
simplified second order energy equation is defined in terms 
of the conductivity function, S = /o(T)dT, which is intro- 
duced so that the energy equation becomes 
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PV 


3h 
z 3z 




l as 

r 3r 



( 1 ) 


where o is the local electrical conductivity, E is the 
voltage gradient, h is the local enthalpy, V z is, axial 
velocity, pV is assumed to be a constant, the flowfield is 
laminar and steady, and the air is treated as being in 
thermodynamic equilibrium. ^ 

As shown in Figure 2a, the stationary arc is assumed 
to take up a large fraction of the constrictor diameter. 

The mass flux is assumed constant, axial pressure gradients 
are assumed negligible and only heat loss due to thermal 
conduction is considered. The analysis takes place only 
within the arc column region, which does not take into 
account flow passing the cathode and anode. The results 
show the axial and radial effects of enthalpy inside the 
column as well as voltage gradients and electrical conduct- 
ivity. Further experimentation has shown that despite the 
shortcomings in the analytical approach, this analysis 
predicts fairly well the energy distribution inside the 
column. 

The first significant analytic thruster model is 
known as the Core-Flow model. The basic difference between 
the Core-Flow theory and the Stine-Watson model is the 
assumption of a thin arc with respect to the constrictor 

\ 
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diameter and the exclusion of axially dependent terms from 
the energy equation. The enthalpy of the core is assumed 
to be given by the Elenbaas -Heller equation and is inde- 
pendent of the mass flow. 5 The basic energy equation, after 
the axial conduction and convection terms are neglected, 
reduces to 


_ F 2 18 dSv 

oE P r r 


where P r is the local radiated power per unit volume. 

As illustrated in Figure 2b, the gas flow in the 
constrictor is broken down into several regions; the hot 
central core where gas is heated directly by the arc, the 
inner flow in which the gas is heated by conduction and an 
outer flow which is heated by contact with the constrictor 
walls. Most of the mass flow occurs in the cool outer flow 


region with negligible flux assumed through the core. 

The poor representation of axial enthalpy flux as 
well as energy flux inside the column at the constrictor 
entrance and exit limit the applicability of the solutions. 
However, the model does give a realistic picture of the 
mass flux at the constrictor inlet and of the overall 
enthalpy distribution.^ The core flow model predicts 
strong gradients in temperature, density and velocity at 
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Figure 2a. Mechanisms of Heat Transfer in 
Stine-Watson Theory ^ 
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performance estimates produced by the Core-Flow model are 
useful approximations and the analysis shows that the 
performance of the arc jet is directly related to the amount 
of flow losses. 

The most extensive numerical solution at that time 

Q 

was done by Watson and Pegot as a comparison to the Stine- 
Watson model. The basic theoretical models were the same 
except fewer simplifying assumptions were made in the 
numerical solution. The analysis of the coaxial flow was 
carried out using a finite difference scheme for varying 
sizes, enthalpy levels and pressure levels. The idea was 
to study the behavior of the wall-constricted arc in the 
range where radiation heat loss is a noticeable portion of 
the total heat loss from the system. Results showed little 
difference between the solutions and that the qualitative 
trends were the same. Watson concluded that the major 
effect of large radiation losses was to lower and flatten 
the enthalpy distribution. 

Another numerical solution of the flow conditions 

9 

inside an arcjet was presented in 1962 by P. Masser , using 
a completely different approach than that of Watson and 
Pegot , and involving the solutions to the equations of 
momentum, energy, continuity, and Ohm's law. The analysis 
assumes the flow is one dimensional, with axial . symmetry 
and that pressure is only a function of x. Axial. heat flow 
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is neglected and all the upstream properties are assumed 
given initially. The results show velocity, enthalpy and 
voltage gradients at various axial locations for hydrogen 
gas and a particular set of input. Recently, Masser used 
an updated version of his program^ to compare the results 

with previous related works, in particular, with Watson's 

8 

numerical results and the experimental results of Todd and 
3 

Sheets . He then studied the effects of varying certain 
parameters, such as velocity and enthalpy. He again makes 
simplifying assumptions to solve the necessary set of 
equations: continuity, momentum, energy and Ohm's law. 

Masser concludes that his solutions are in reasonable 
agreement with Todd and Sheets and shows that his voltage 
and enthalpy gradient profiles closely match those of 
Watson. Results also indicate an increase in performance 
of the engine and better correlation with other results as 
the arc chamber length is increased. 

1.3 Current Analysis 

A new approach to solving the two dimensional 
Navier-Stokes equations numerically has been developed 
recently by C. Merkle, G. Molvik, and D. Choi^. The 
original application was to solve a compressible or incom- 
pressible laser heat addition problem. For the case of 
the arcjet, the geometry was adapted to represent- the 


10 


constrictor section of the arcjet and the heat addition 
modified to represent a concentrated source on the center- 
line of the constrictor. Analytical and numerical work on 
this problem to date has not completely included the effects 
of viscosity in the equations. The code can be modified to 
model several different boundary conditions that depend on 
which specific initial conditions are being used. 

The solution consists of formulating the equations 
into matrix form and applying an Euler-Implicit finite 
difference scheme. The gas flow inside the constrictor is 
taken as laminar and only steady state solutions are con- 
sidered. A power source term is added to the energy equa- 
tion for the arc, viscosity is taken into account, and real 
gas properties are used in solving for the local thermo- 
dynamic properties inside the constrictor. The energy 
input is assumed constant in the axial direction, falling 
off radially in a gaussian distribution to zero at some 
specified distance from the centerline. The heat is 
as siimed to start at the tip of the cathode and stop at the 
constrictor exit. The region covered in the analysis and 
the representative finite difference grid spacings are 
shown in Figure 3. The steep gradient of cool to hot gas 
in the constrictor produces a highly two dimensional flow- 
field. 12 




Chapter 2 


FORMULATION OF THE PROBLEM 


2.1 Governing Equations 


The two-dimensional unsteady compressible Navier 

12 

Stokes equations are used to describe the flow in the 
constrictor. In cartesian coordinates the conservation 
vector form of these equations is 


3^ , 3E , 3F 
3t 3x 3y 


3V 

3x 


3W 

3y 


+ H 


(3) 


where 



(PU 

2 

pu + p 
puv 

i(e + p)uy 


F = 




In addition, p = density, u = axial velocity, v = radial 
velocity, p * pressure, e = internal energy, and q = heat. 

In this analysis, p, pu, pv, and e are the dependent 
variables, which are contained in the vector Q. The equa- 
tions are solved in terms of these parameters, namely 
continuity, x-momentum, y-momentum, and energy. The flux 
vectors E and F represent the convective terms, V and W are 
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the viscous terms and H is the source term. The viscous 
terms incorporate all shear stresses and there is a pres- 
sure dependence in both the x and y directions implying 

12 

full consideration of. the pressure gradients . An ex- 
panded component form of the above governing equations is 
contained in the appendix. 

Pressure is obtained from the equation of state 

where p * density, R « universal gas constant, T = absolute 

temperature and = molecular weight. 

The source term added to the energy equation is a 

2 

guassian distribution, (P = ^), that is specified as 

taking up a certain percentage of the constrictor diameter. 
The percentage of the total diameter that is covered by the 
arc is defined as Percent. The quantity Py is subsequently 
defined from this as 


Py = (1 = Percent )Y, 


max 


( 5 ) 


where Y is the constrictor radius. Ther term b is 
max 

defined as 


log(i 


b = 


max 


( 6 ) 


(Py - Y ) 
' 3 max 7 
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and finally, the power P becomes 

p «= p e^^max 7 ^ 

max 


(7) 


where P max is the specified maximum power input along the 
centerline, y being the radial coordinate above the wall of 
the constrictor. This is empirically shown in Figure 4. 

The power is added directly to the energy equation in SI 

3 

units of W/m , where the third dimension is considered to 
be of unit length and the energy input remains constant 
through each iteration. 

The equations are written in a generalized coordi- 
nate system in order to simplify the curvilinear grid gener- 
ation. Although this is not necessary for the particular 
case studied here, any future adaptations to the shape of 
the constrictor can be obtained with only minor modifica- 
tions to the geometry specifications in the code. 

2.2 Euler-Implicit Method of Solution 

The governing equations are solved using a time 
dependent finite difference scheme. The term Euler-Implicit 
means that the first order differences with respect to time 
are implicit and the solutions are unconditionally stable. 
The partial derivatives with respect to space are then 
second order differences. Viscosity and the power source 

are handled explicitly throughout in the terms on the 
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right hand side of the equations. All terms on the left 

hand side of the equations are solved implicitly. The 

equations are solved in conservation form so that mass , 

momentum and energy are conserved identically and linear- 

13 

ization is required only with respect to time. 

As stated previously, the vector form of the Navier 
Stokes equations is 


3^ . 3E , 3F 
3t 3x 3y 


3V 

3x 


3W , 

3y 


H 


( 8 ) 


Solving these equations simultaneously using the 
Euler-Implicit scheme requires writing them in delta nota- 
tion. To do this, we first define the Jacobian matrixes A 
and B as 



and 



The original delta form of the equation is 


(9) 


(I + At + At -|^)AQ = -At[-^(E-V) + -^(F-W) - H] (10) 

Since the terms on the right hand side of the equa- 
tion are solved explicitly, they can be referred to simply 
as the quanity 

\ 
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RHS * -At[-^(E-V) + ~(F-W) - H] (11) 

The Douglas-Gunn approximate factorization^ is 
applied to the equations at each time step. The implicit 
terms are written in the form 


(i + At|£)(I + At|5)AQ - RHS (12) 

AQ* is now defined as 

AQ* = (I + At||)AQ (13) 

and then 

AQ* = (I + At||) _1 RHS (1A) 

so the above equation can be rewritten 

(I + At||)AQ* = RHS (15) 

The desired quantity AQ can now be written as 

AQ - (I + ||)’' 1 AQ* (16) 


Blottner’s block solver"^ is used to solve the block 
tridiagonal matrix first for the vector AQ* and then once 
again for AQ. 
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The quantity AQ is used as the basis for computing 
the values at the next time step. The quantities making up 
AQ are Ap, Apu, Apv, and Ae. These quantities are added to 
the previous value at that location and all thermodynamic 
properties are computed from these relationships. 

A constant time step is used in the iterations for 
simplicity. A variable time step produces a faster conver- 
gence, but is not necessary for. the steady state solution 
used. The time step size is a function of the CFL condi- 
tion for stability and is limited by the value of the CFL. 

If the CFL is too large or too small, the solution will not 
converge. The time step is obtained from the inputed CFL 
and the maximum eigenvalue from the Jacobian matrices A and 
B, which is a function of the speed of sound. The CFL value 
used is 5, which is considered optimum by experimentation. 

2.3 Boundary and Initial Conditions 

The most important criteria for correctly solving 
the Navier-Stokes equations using the Euler- Implicit scheme 
is to have valid boundary conditions. The finite differ- 
ence grid for the constrictor diameter is generated first in 
rectangular coordinates from the specified number of x and 
y points needed for an accurate analysis, and then trans- 
formed to generalized coordinates £(I,J), The J=1 

location corresponds to the lower constrictor wall, J=JL 

j 
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represents the upper wall, 1*1 is the inlet and I=IL is 
then the constrictor exit, as Figure 5 shows. 

The initial upstream conditions specified are stag- 
nation pressure, stagnation temperature, and a variable 
velocity distribution. Downstream, a constant pressure is 
imposed at the constrictor exit. The code was determined 
to work best for cases with a constant velocity profile and 
slip conditions at the wall, which corresponds to higher 
Mach number, compressible flow conditions. The constant 
velocity profile is defined as u*constant along the entire 
inlet location. No radial component in the initial 
velocity is assumed, which is written in the form arctan 
v/u = 0. 





Chapter III 


NUMERICAL RESULTS AND DISCUSSION 
3.1 Cold Flow Results 

Convergence is most important in an iterative 
solution over time. The solution is not completely correct 
until the dependent variables no longer change with respect 
to time. These parameters, Ap/p, Apu/pu, Apv/pv, Ae/e, 
become zero when the solution has reached complete conver- 
gence. The system is then considered to continue in steady 
state. A solution diverges when one or more of the input 
data are incompatible with the set boundary conditions of 
the problem. The input conditions in this analysis that 
can be varied to simulate a certain system are the inlet 
stagnation temperature (T 0 ) , the inlet stagnation pressure 
(P Q ) , the back pressure, the amount of energy input per 
cubic meter of the arc (Power) , and the physical size 
of the constrictor section. In addition, significant 
differences in the solution occur with changes in the 
grid size. The number of iterations required to reach 
a converged solution changes with the grid size and to 
a lesser extent with the specific initial conditions. 
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Simple test cases of the program were run to deter- 
mine if a converged solution was obtainable with the modi- 
fied heat input and initial conditions. A solution is 
easily predicted for a uniform flow with no heat addition; 
with a constant velocity profile, slip conditions at the 
walls, and uniform initial temperature and pressure 
radially, the flow should remain unchanged because it is 
already in steady state. After several different sets of 
initial conditions were run, the flow characteristics 
remained constant with time and all convergence parameters 
began and remained at zero. These results were consistent 
throughout any variation in geometry or grid spacing. 

The final range of initial conditions used in this 
study are similar to those successfully used in other 
applications of this code and are representative of real 
situations that might occur inside a constricted arc. 
Stagnation temperatures range from 800° to 3200°K, inlet 
Mach numbers run between .18 and .5, and stagnation pres- 
sures are between 1 and 6 atms. The physical size of the 
constrictor sections modeled ranged between lengths of .1 
and .25 m and heights between .025 and .1 m, however, 
throughout the study it was noticed that the actual physi- 
cal size of the constrictor channel did not affect conver- 
gence. The number of grid points chosen in the x-direction 
is in every case much smaller than the number of ^points 
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chosen in the y-direction because the property profiles 
change much more drastically across the diameter than from 
one axial location to another. Faster convergence occurs 
with a coarser grid mesh, however in most cases a finer 
grid mesh was used because of the increase in accuracy of 
the solutions. 

The parameter with the greatest effect on the con- 
vergence of the solution is the power input. The actual 
percentage of area covered by the arc and its corresponding 
percentage of the maximum possible power in that area 

determine how many Watts of power are emitted into the 

2 

constrictor. The power inputs run between 1.2x10 W and 
1.8x10^ W, or 1x10^ and 1.5x10^ W/cm^, where the value 
stated represents the level of power at the centerline 
which falls off in a gaussian distribution to zero at a 
distance of 20 percent in the constrictor diameter on 
either side of the centerline. 

3.2 Effect on Power Input 

Once it was successfully determined that the code 
was valid for a uniform flow with no heat source, energy 
could be added in small increments to determine the effect 
of the presence of the arc on the flow characteristics of 
the moving gas and convergence of the solution. The solu- 
tions to the code are outputted in the form of the values 
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of the local thermodynamic properties at each point on the 
grid after the desired number of iterations. The most 
important of these parameters is temperature, which is a 
direct indication of how much energy is being absorbed by 
the hydrogen gas from the arc. 

The lowest power input to record a significant 

a 2 

temperature rise was 1x10 W/cm where the initial condi- 
tions used were a stagnation temperature T=1500°K, a stagna- 
tion pressure of P=3 atm, and an inlet Mach number of 
M=.254, using a grid size of 11 by 20. The actual size of 
the constrictor was .15 m long by .05 m high. The solution 
quickly converges to zero and the maximum temperature rise 
was 2°K, which occurs on the centerline at the exit. Con- 
vergence was also computed for power levels of 1x10^, 

1x10^ and lxlO^^W/cm^. Figure 6 shows a log-log plot of 
power input vs. AT and indicates that as the power increases 
by a factor of ten, the maximum temperature rise also in- 
creases by a factor of ten. Convergence plots for these 
cases indicate that the lower the power input, the faster 
the convergence. This implies that the lower the heat 
addition, the less disruption is introduced into the systems 
and the faster the system can reach a steady state. 

The converged solutions were checked by determining 
if energy was conserved in the system. The energy at the 
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constrictor should equal the energy at the exit. The 
program has internal limitations as to the accuracy of the 
comer grid points at the inlet and exit and the calcula- 
tions were performed with increased accuracy at axial loca- 
tions one point after the inlet and one point before the 
exit, 1-2 and I-IL-1. 

The energy difference between the inlet and the exit 
is calculated by solving the following equation, 

Q “ / (e + p)udy - f (e+p)udy (17) 

out in 

where Q is the power of the arc in Watts . A program was 
designed to take the converged solution variables, e, p, and 
u and automatically integrate the inlet and exit energies . 
The quantity Q is determined by the amount of area in the 
entire constrictor that is actually covered by the arc. 

These calculations show consistent error of 6 to 7 percent 
between the difference in power from the inlet. to exit and 
the power inputed from the arc. This is considered reason- 
able because errors may be introduced in the exact calcula- 
tion of Q. 

Although there is no lower bound heat input for this 
code, there is an upper limit to the amount of energy that 
can be introduced into the system and still reach a con- 
verged solution. For a given set of initial conditions, 

i 
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when the power per square meter is raised past a certain 
value, the dependent parameters, AQ/Q, become so large that 
they caused an arithmetic fault in the program. The limit 
was found to be somewhat dependent on the fineness of the 
grid spacings but overall the power level could not be 
raised higher than about 1x10^ ^W/cm^ corresponding to a AT 
of about 2000*R. 

Figure 7 shows the temperature cross section at the 
axial location just before the exit, where the grid points 
were 11 x 40. For given conditions of initial stagnation 
temperature T=2000°K, stagnation pressure P=4 atm, constant 
inlet Mach number M=.19 and Power=lxl0^ 2 W/ cm 2 , the maximum 
AT is 1800 degrees. The cross section is similar to the 
gaussian power input profile (Figure 4) except appears 
wider. As expected, the temperature is highest where the 
concentration of energy is greatest. Figure 8 shows the 
log-linear plot of the convergence criteria for this case. 
Ap/p, Apu/pu, Apv/pv and Ae/e become zero as the steady 
solution approaches . The AQ/Q terms do not smoothly 
approach zero, but rather take short jumps up and down 
while gradually decreasing toward zero. In every case run, 
the x-momentum took the most iterations to converge to zero. 
The program was run again, increasing the power input to 
1.3xl0 12 W/cm 2 and after 165 iterations, the AQ/.Q diverged 
and the code automatically stopped. The code can 'allow only 
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so much energy input into the constrictor before the energy 
equation registers a discontinuity which affects the energy 
convergence and ultimately the rest of the convergence 
parameters. However, it is interesting to note that the 
solution begins with AQ/Q decreasing normally, then suddenly 
within only a few iterations the corrections increase 
dramatically and the solution fails . Figure 9 shows the 
drastic change in AQ/Q. The values jump to hundreds of 
times their initial values before the code stops because of 
an arithmetic fault. 

3.3 Effect of Temperature Initial Conditions 

Initial inlet and wall temperatures have an effect 
on the maximum overall change in temperature of the con- 
strictor for a given value of energy input. Six cases were 
run with all the same conditions except the initial temper- 
ature was varied. Figure 10 shows the temperature profiles 
at the exit axial location. The initial conditions used 
are P=4 atm, and T=400°, 800°, 1400°, 2000°, 2600°, and 
3200°K. Note that as the inlet temperature increases, the 
maximum change in temperature on the centerline decreases, 
but in each case the profiles have the same shape. When 
T=800°K, the AT was 2300 degrees whereas with 3200°K, the 
maximum temperature change was only approximately 800°K. 

The solution diverged, however, when run with temperatures 

j 












Y-Coordinate 


Figure 11. Temperature Cross Sections at Various Axial Locations 
(T=2000°K, P=4 atm. M=.35, Power=1.5X10l2W/ cm 2) 
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below 400°K, implying that there is at least a lower limit 

for temperature in the range of obtainable solutions. 

Figure 11 shows the axial change in temperature 

profiles for the case of T=800°, P=4 atm, and M=.348, with a 

12 2 

power input of 1.5x10 W/cm . The local temperature cross 
section is plotted for x locations of IL=1,3,5,7,9 and 11, 

where IL-1 corresponds to the inlet and IL-11 is the exit. 

The maximum temperature for this completely converged 
solution occurs at the exit, as well as the maximum 
velocity. 


3.4 Effect of Velocity and Pressure 

Both velocity and pressure changes affect the conver- 
gence and temperature rise in the constrictor. The velocity 
distribution does not remain constant after the heat source 
is introduced into the system. Figure 12 shows the velocity 
distributions for the same initial conditions and axial 
locations as the above temperature profiles in Figure 11. 

The velocity starts out as a constant across the inlet but 
as the program is run and heat is added, the velocity 
decreases in the center region due to interaction with the 
arc. The velocity increases fairly steadily as it moves 
downstream until right before the exit where it takes an 
unprecedented jump at the center while forming wells of low 
velocity points on either side of the centerline. At the 
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Figure 12. Cross Sectional Velocity Profiles at 
Various Axial Locations 
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11 2 

same time, with heat inputs above 1x10 W/cm , the 
velocity gradients in the y-direction become significant, 
causing the flow to become very turbulent in the center 
region. However, the velocity remains radially uniform 
throughout the channel outside the region of the arc. 

Pressure follows a similar but opposite pattern to 
the velocity distribution as can be seen in Figure 13. 

The pressure begins by rising at the inlet in the center of 
the constrictor but decreases consistently and forms a more 
linear profile as it moves axially downstream, decreasing 
over 10 percent from its initial value of 4 atm when it 
reaches the exit . 

A lower limit is also apparent for possible initial 

stagnation pressures for a specific set of input. For 

example, one case was run with an initial temperature 

T=1500°K, an inlet Mach number M=.254, and a power input 
11 2 

Power=5xl0 W/cm , where the physical size of the cons trie 
tor was .1 m long by .025 m high. The solution with an 
initial stagnation pressure of 1 atm diverges, while the 
solutions with pressures of 6 and 12 atms converge per- 
fectly. The same holds true for a lower limit of velocity 
where a velocity of 400 m/sec diverges and 600 m/sec gives 
a converged solution. 
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Chapter IV 


CONCLUSIONS AND RECOMMENDATIONS 

This paper has described how the two dimensional 
compressible Navier -Stoke' s equations can be solved using 
an Euler Implicit finite difference solution, with heat 
addition and boundary conditions representing the constric- 
tor section of an arcjet. The code produces converged 
solutions for varying boundary and initial conditions , and 
can show how changing one or more of these conditions can 
effect the solution. 

The temperature profiles prove that the arc energy 
source does indeed heat up the hydrogen gas passing through 
the arc. The maximum temperature rise occurs along the 
centerline where the concentration of energy is the greatest. 
The overall temperature profile throughout the constrictor 
closely resembles the gaussian shape of the arc heat inten- 
sity. As the heat input is increased, the maximum change 
in temperature from the initial conditions also increases. 
The velocity of the gas flow along the centerline decreases 
at the inlet as a result of interaction with the arc but 
rises steadily until the exit where the profile suddenly 
becomes highly inflected. However, outside of the arc 
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column region, the velocity profile remains uniform and 
increases linearly. The corresponding pressure effects are 
opposite except the profile maintains relative uniformity 
radially while decreasing in the axial direction. 

Although the code produces theoretically very 
accurate results because few simplifying assumptions are 
made in the solution to the original equations , the range 
of initial conditions that give converged solutions is 
prohibitive for extensive research in the area of con- 
stricted arcs with large power sources. The fact that the 
program cannot give a converged solution for power inputs 
greater than 1x10 W/cm is a major limitation to the code 
but it is likely that further research can solve this 
inadequacy and also provide a code that can produce reliable 
results in the incompressible velocity range. The Euler 
Implicit approach to numerically solving the arc jet problem 
is valid, however more work is needed to extend the range 
of conditions that will produce converged solutions. 
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APPENDIX 


Expanded Form of the Governing Equations 


Beginning with the two dimensional vector form of the 
governing equations 


3^ , 3E , 3F 

at ax ay 


av . aw . „ 
a£ + . ay + H 


Expanding each of the vectors Q, E, F, V, and W gives the 
four equations of continuity, x-momentum, y-momentum and 
energy . 
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3p . 

at 

a ( pu) . a(pv) 
ax ay 

= 0 


a (pu) 
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XflH + iVv , ( 3u + 3 Vn + 2u— 

M 3x 3y J M 3x 3y' y 8y 


+ + + 


+ x <i + + 2 u( li u + + K I + if> 


where K = thermal conductivity, y = viscosity, and 
X = second coefficient of viscosity which are tabulated 
values . 


42 



End of Document 



